The Annals of Applied Probability 

2006, Vol. 16, No. 1, 107-139 

DOI: 10.1214/105051605000000656 

© Institute of Mathematical Statistics, 2006 

A SCHEME FOR SIMULATING ONE-DIMENSIONAL DIFFUSION 
PROCESSES WITH DISCONTINUOUS COEFFICIENTS 

By Antoine Lejay^ and Miguel Martinez 
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The aim of this article is to provide a scheme for simulating difltu- 
sion processes evolving in one-dimensional discontinuous media. This 
scheme does not rely on smoothing the coefficients that appear in the 
infinitesimal generator of the diffusion processes, but uses instead an 
exact description of the behavior of their trajectories when they reach 
the points of discontinuity. This description is supplied with the local 
comparison of the trajectories of the diffusion processes with those 
of a skew Brownian motion. 



1. Introduction. The aim of this article is to provide a scheme for the 
simulation of the one-dimensional diffusion process X generated by the dif- 
ferential operator 

(1) L = -—[a—\+h—. 

^ ' 2dx \ dx) dx 

Here, a, p and b denote piecewise smooth functions that may have an infi- 
nite number of discontinuities on a countable set of points J . We assume, 
however, that J has no cluster points and that the functions a, p and b 
have everywhere left and right limits. The triplet (a, p, b) will be called the 
characteristic of L. 

If a belongs to C^, L may be transformed into 

ap d'^ f a' \ d 

and we see that, even if the coefficients are smooth, using a Euler scheme for 
the simulation of X requires us to compute the derivative of a, which may 
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be quite expensive from a numerical point of view. However, if a = 1, it has 
been proved in [39] that the Euler scheme converges but, yet, specifying its 
speed of convergence still remains an open and challenging problem. 

The scheme presented in this article is different from the Euler scheme. 
It should rather be seen as a variation of the well-known random walk on 
spheres. 

The basic idea is the following. First, we replace the differential operator L 
by another one whose coefficients are piecewise constant, which provides 
good approximations of the solutions of the elliptic and parabolic PDEs in- 
volving L (see Section 8 for a computation of the error). Second, we simulate 
the stochastic process generated by the approximation of L at a given time 
using a description of its behavior when it is around a point where a or p (or 
both) are discontinuous. Mainly, we compute quantities related to the first 
exit time and position on some intervals. This method is exact in the case 
of piecewise constant coefficients when 6 = because it describes correctly 
the behavior of the diffusion process when it reaches a point in J . 

Let us explain our approach with a simple example: assume that 6 = and 
suppose that {a{x),p{x)) = (a+, p"*") if x > and that {a{x),p{x)) = (a~ , p~) 
if x < 0. Here, , are positive constants. Let u be the weak solution of 
the parabolic problem associated with L: 

(2) ^t^^hA = Lu{t,x) and u{Q,x) = f{x). 

It is well known that this problem is equivalent to the following transmission 
problem (see [18], e.g.): 

1 ^!f^ = ^i^M^)An(t,x), onM; xM; andonM; xMl, 
[a+Vu(i,0+) =a-Vn(t,0-), for t > 0. 
Now, let us introduce ^ defined as follows: 

X X 
^(^) = r^^ l{x>0} + , l{x<0}- 

Since the function u is of class on ]R_|_ x M \ {0}, it is easy to check that 
v{t,x) := u{t, ^{x)) is the solution of another transmission problem 

( dv{t,x) ^ 1^^^^^^^)^ on X and on X Ml, 



" =Vt;(t,0+) = — ^Vi)(t,0-), fori>0. 



VP'^ VP" 
(4) 

Formally, this new transmission problem (4) is also equivalent to the parabolic 
problem ^ = Lv{t,x), where L is the formal differentiable operator + 
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(55qV, with 



Va(0+)/p(0+) - ^a(0-)/p(0)- 
Va(0+)/p(0+) + ^a(0-)/p(0-) 



As shown in [30, 31] (see also [11]), this differentiable operator L is the 
infinitesimal generator of the skew Brownian motion. This process can be 
constructed from a reflected Brownian motion by simply choosing the sign 
of each excursion by tossing an independent Bernoulli random variable of 
parameter {(3 + l)/2. Heuristically, this means that the particle chooses to 
go on (resp. M_) with probability {(3 + l)/2 [resp. (1 — /3)/2] when it 
reaches 0. Unfortunately, this description is not relevant due to the fact that 
there are infinitely many small excursions. Besides, this approach does not 
permit us to understand the real behavior of the particle as a function of 
and . However, numerical simulation of the skew Brownian motion is easy, 
and using a simple and deterministic change of scale, it is possible to solve 
(2) using the probabilistic representation u{t,x) =Kx[f{^~^{Z^))], where 
is a skew Brownian motion of parameter /3. 
Instead of working with PDEs, one may prefer to use the Ito-Tanaka 
formula as in [29]. This can be achieved using a precise description of the 
process X related to L. As shown in [20], Chapter 5, X solves the following 



where i? is a Brownian motion and is the symmetric local time of X 
at 0. The diffusion process ^{X) is then a skew Brownian motion with the 
coefficient /? given above. 

The basic idea of this paper is to use this kind of description because of 
the particular properties of the skew Brownian motion. 

Although this is not the first use of the skew Brownian motion in Monte 
Carlo methods (see [21, 40] or, more recently, for financial applications, [9]) 
or in modeling (see [7] for application in ecology), there has been neither a 
systematic study of this process in this framework, nor a complete exposition 
of the interplay between the different coefficients. 

The numerical method presented in this article follows the idea of [22], 
where diffusions with constant coefficients on the edges of a graph were sim- 
ulated. As the infinitesimal generator provides locally the behavior of the 
particle, some properties of the skew Brownian motion allow, in this partic- 
ular framework, to describe what happened at the nodes of the graph. These 
are the main reasons which explain our choice of approaching nonconstant 
coefficients using piecewise constant approximations rather than smoothing 
the discontinuities around the points of J^. Indeed, this last procedure turns 
out to be unstable in practice and very expensive numerically. 



SDE: 
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We will also provide results concerning the diffusion process generated 
by L under more general assumptions than in [20], Chapter 5, mainly by 
describing it as the solution of some SDE involving the local time and being 
aware of the boundary conditions. In his Ph.D. thesis [27], one of the au- 
thors gives other results on diffusion processes generated by divergence-form 
operators, and uses also space and time transforms as a very natural tool. 
Among his results, he studies the speed of convergence of the Euler scheme 
with discontinuous coefficients obtained after making use of some other scale 
transform. 

One could also think of using a random walk as proposed in [23] (after 
a proper change of the scale). The advantage of this method is that the 
time step is incremented with a constant value and not with a random 
variable. This perspective is studied by P. Etore in the recent article [12]. 
We also argue that our algorithm can be implemented locally around the 
points where discontinuities hold. One may use more efficient algorithms 
(Euler scheme, Milstein scheme, . . . ) in the regions where the coefficients 
are smooth. 

Outline of the article. In Sections 2 and 3 we show how to construct 
the stochastic process generated by L. In Sections 4 and 5 we study the 
properties of the semi-group generated by L and the convergence of the 
solutions of the PDEs with respect to a family of differential operators. In 
Section 6 and 7 we study the process X generated by L and we show how to 
transform it into a process that behaves locally like a skew Brownian motion. 
The algorithm is explicitly given in Section 9 and its error is studied in 
Section 8. Finally, as an example, we show numerical results concerning the 
density of a doubly skewed Brownian motion and we compare this scheme 
with others for a differential operator with nonconstant coefficients. 

Hypotheses. Let £ < r be two numbers that belong to M. It is possible 
that H. = —oo and/or r = +oo. 

Assume that the coefficients a, p and b are defined on [i,r] and satisfy 

(5a) a, p and b are measurable. 



for some constants A, A > 0. 

Let ^{dx) be the measure p{x)~^ dx on (•£,r) and G denote the open set 
(£, r). It is possible that G = M. For technical reasons, we restrict ourselves 
to the case where ^ and r are simultaneously finite or infinite. In fact, all 
the results given here can be extended for G = {I, +oo) or G = (oo, r). 



(5b) 
(5c) 



for all X G [i,r], p{x) E [A, A] and a{x) G [A, A] 
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We define L^(G) as the space of measurable functions on G that are 
square integrable with respect to the Lebesgue measure. We also define 
H^(G) [resp. Hq(G)] as the completion of the space of smooth functions 
(resp. smooth functions with compact support) on G with respect to the 
norm 



HHa) = v/|l/llL^(G) + llV/ll^.(a)- 



Recall that, for any connected interval G, all functions in the Sobolev 
space H-^(G) have a continuous version. In all the following we will system- 
atically identify a function in H^(G) with its continuous version. 

2. Removing the drift. We assume that a, p and b are smooth. 
When G is finite, choose (L,Dom(L)) to be any of the differential opera- 
tors: 

J L is given by (1), 

1 Dom(L) = {/ G C2(G;M)|/(r) = /(^) = 0} 
for Dirichlet boundary condition (b.c.) or 
J L is given by (1), 

1 Dom(L) = {/ G C2(G;M)|/'(r) = /'(£) = 0} 

for Neumann b.c. 

When G = M, let (L,Dom(L)) be the differential operator: 

J L is given by (1), 
\Dom(L) =C^(M;M). 

Because a, p and h are assumed to be smooth, it is well known that (L, Dom(L)) 
is the infinitesimal generator of a continuous strong Markov process {X, (Px)x6g) 
(see [5], e.g.) and, in addition, this process is the solution to the stochastic 
differential equation (SDE) 

dXt = ^{Xt) dBt + + ^) i^t) dt. 
If 6 = 0, then (L,Dom(L)) can be associated to a symmetric bilinear form 
£{u,v) = ^ J a{x)u {x)v' {x) dx 
defined for u,v in Cj!,(G;M) through 
(6) £iu,v)=lLu{xHx)p~\x)dx 



for (n,z;) GDom(L) x Cj^(G;M). 

Since the symmetry of L will play an important role in most of the re- 
sults given in this article, we would like to be able to transform (2) into an 
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equivalent form without drift: as we will now see, this is possible thanks to 
the crucial fact that we are working in the one-dimensional space. Let us 
explain one way of removing the drift by transforming a and p: let ^ be the 
function 

(7) ^{x)=rh{x)dx with h{x) =2 T \ , dy. 

Jo Jo p{y)a{y) 

If ^' is given by (7) and is bounded, a simple calculation shows that 
2 dx\ dx J 2 dx\ dx ) dx 

and we see that 

/ON -^P d ( ^ d\ 

(8) e -—[ ae — ] 
^ ^ 2dx\ dx) 

is a linear operator which is symmetric with respect to the measure p{x)~^e^ dx 
and, thus, we can manipulate (6). 

If G is bounded, then ^ is bounded by a constant that depends only on A, 
A and the Lebesgue measure Meas(G) of G. Hence, this procedure is valid. 

If G is not bounded, then ^ is not bounded in general, and we cannot 
make this transformation without using further arguments. However, we 
claim that it is possible to replace ^, at least locally, with ^ defined, for 
example, by 

= ^(x) - ^{n) if x G [n, n + 1] for all n G Z. 

Thus, it is possible to repeat our argument and to transform L locally in a 
linear operator which is symmetric with respect to the measure p{x)~^e^^^^ dx. 

As a matter of fact, instead of considering {a,p,b) as the characteris- 
tic of L, one should only consider L related with the new characteristic 
(e*a, e~*p, 0). 

Note that the transform used here to remove the drift applies for nons- 
mooth coefficients a, p and b since it is always possible to use a regularization 
procedure and pass to the limit: see regularization results in Section 5. 

3. Existence of a stochastic process. In dimension one, results on the 
existence of stochastic processes generated by (L,Dom(L)) may be proved 
using (at least) three ways: using the properties of the density transition 
function, using the scale function and the speed measure, or using the Dirich- 
let form theory. We will see that all these methods lead to the same process. 
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3.1. Using Dirichlet forms. Let us assume at first tliat 6 = 0. This will 
allow us to use the theory of symmetric Dirichlet forms. We have seen in 
Section 2 that this hypothesis is not a restriction. 

Let £ be the bilinear form 



^ 2 Jg dx dx ' 
on 'L?{G;i'{dx)). The domain Dom(i5) is 

Dom(iS) = H^([£, r]; ;/((ix)) ~H"'^([^,r]) for the Neumann b.c, 



It is well known that {£,Y)om.{£)) is a regular, local Dirichlet form [13]. 
Hence, it generates a continuous strong Markov process [X, {J-'t)t>o, i^x)xeG)- 
Besides, the process is conservative if G = M or Dom(iS) = H^([^,r]) which 
corresponds to the Neumann b.c. See Lemma 1.6.5 of [13]: recurrence fol- 
lows by applying Theorem 1.6.3 of [13] and the fact that, in both cases, 
1 e Dom(£:) and £{1, 1) = 0. 

Besides, if G = {£, r) and Dom(£:) = Hj(G), the Dirichlet form {£, Dom(£:)) 
still possesses the local property: thus, we can repeat all the arguments of 
Example 4.5.1, page 166 in [13] and X is absorbing at both end-points 
of G (see also Theorem 4.2.2, page 154 in [13]). We are allowed to write 
Xt = 1|t-<j}.5 -|- l|T->(}.Yj, where 6 is the "point at infinity" added to M, 
Y is the process generated by (£^,Ho(M)) with a = p = l outside G, and 



Remark 1. As the dimension of the space is one, each point x G [i,r] 
is nonpolar and has a positive capacity, and all the statements of type "for 
quasi-every point of [i,r]" mean in fact "for every point of [^,r]." 

Remark 2. One could assume that the coefficients a and p are locally 
bounded and locally uniformly elliptic. In this case, the process is not nec- 
essarily conservative, which means that it may explode in finite time. This 
issue is discussed in [35]. 

3.2. Using the properties of the semi-group. To the Dirichlet form, {£, 
Dom{£)) may be associated a linear operator (L,Dom(L)) on L'^{G;i'{dx)) 
defined through the relation 

£{u,v) = —{Lu,v)i^2(Q)v(dx) foi' i'^iv) € Dom(L) x Dom(<f). 
The operator L is the one given by (1) with domain 




for all u,v & Dom{£) 



Bom{£) = RliG; u{dx)) ~ RI{G) 



for the Dirichlet b.c. 



T = mi{t>0\Yt(^G}. 



Dom(L) = {/ G Hj(G)|L/ G L'^{G)} for the Dirichlet b.c, 

Dom(L) = {/ G H^([^,r])|L/ G L^{G)} for the Neumann b.c. 
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It is possible to show that the operator (L,Dom(L)) is the infinitesimal 
generator of a semi-group {Pt)t>o- We will see in Section 4 that Pt has a 
density transition function p{t,x,y) and that some general estimates hold 
for it. These estimates [see (12), e.g.] are sufficient to ensure the existence 
of a strong Feller, continuous process {X, {J-'t)t>o, {^x)xeG)- 

3.3. Using the scale function and the speed measure. Using the scale 
function and the speed measure gives another way to define the stochastic 
process {X,{Tt)t>o,iPx)xeG)- 

Let us set, for x in {i,r), 

(9) h{x) = 2 r ^^[y] , dy, m(dx) = ^^^^fix, 
^ ^ ^ ^ Jo p{y)a{y) ^ ' p{x) 

(10) S{x) = r ^^P(~^(^)) and V{x)= r m{x)dx. 

Jo a[y) Jo 

For the Dirichlet b.c. (strictly speaking, the choice of m in this case is that 
of an absorbing condition), we set m{{y}) = +oo for y = £ or y = r. For 
the Neumann b.c, m{{y}) = for y = r or y = £. On that topic, see, for 
example, [5] or [34]. We will see in Corollary 1 that this process corresponds 
to the one already constructed using Dirichlet forms. 

Unless a is smooth, the process X is not a semi-martingale in general. 
It is a Dirichlet process. Nevertheless, some stochastic calculus for X is 
possible using the theory of Dirichlet forms and time reversal techniques: 
see [13, 26, 32], for example. 

4. Properties of the semi-group. Let G be the open set G = {£, r) or 
G = M. For a function u in H^(G), we set Tu{x) = u{x) if we are interested 
in the Dirichlet b.c. and Tu(x) = if we are interested in the Neumann 

b.c. Note that ^^^^ is a distribution in general and might not be well defined 
at all points, but in this section through abuse of notation, we will use this 
symbol as a notation for the distributional derivative of u. 

Let us consider the parabolic partial differential equation (PDE) 



' du(t,x) p(x) d ( . .du(t,x) 
— a[x) 



(11) 



dt 2 dx\ dx 

+ on(0,oo)xG, 

Tn(t,x)=0, on (0,oo) X {^,r}, 

,u{0,x)=ip{x), xGG. 



Let us consider the parabolic PDE (11). Unless a, p and b are sufficiently 
smooth, the solution u is a weak solution that can only be chosen in the 
space C(0,r;L2(G)) n L2(0, T; HJ(G)) in the case of Dirichlet b.c. and in 
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C(0,r;L2(G))nL2(0,r;Hi([^,r])) in the case of Neumann b.c. Such a solu- 
tion exists and is unique as long as belongs to L^(G). 

It is standard that (L,Dom(L)) is the infinitesimal generator of a semi- 
group (Pt)t>o on L'^{G,iy). 

Proposition 1. (i) For any t > 0, Pt has a positive density function 
p{t,x,y) with respect to the measure v. Besides, 

u{t,x)= / p{t,x,y)ip{y)diy{y) 



G 

is a version of the solution of (11) which is continuous with respect to {t,x) 
on (0, oo) X G. 

(ii) The function {t,x,y) p{t,x,y) is {a/2,a,a)-Hdlder continuous in 
{t,x,y) on every compact of (0, oo) x G, where a depends only on X, A and 
the size of the chosen compact. 

Proof. The existence of a density of Pt is a classical result in the theory 
of PDEs. A proof of (ii) may be found, for example, in [2] . □ 

Proposition 2. (i) If G = R or G is bounded and the Dirichlet b.c. are 
used, then there exists constants Ci and C2 depending only on X, A and T 
such that 

(12) p{t, X, y) < -j== exp 



/27rt V t 

for any {t,x,y) G [0, T] x M x M. If b = 0, then the estimate (12) is uni- 
form in T. This bound is called Aronson^s estimate. [Note that, however, D. 
Aronson proved these estimates in a general case, but it was initially proved 
simultaneously for operators of type V{aSJ-) by E. De Giorgi and J. Nash.] 
(ii) // G is bounded and the Neumann b.c. are used, then, for any T > 
and any 9 > 1/2, there exists some constants Gi, G2 depending only on X, 
A, 9, r, i and T such that 

(13) p(t,a;,y) < ^exp'^'"^^'^ 



t^ t 
for all t G (0, T] and all x,y [i, r] . 

Proof, (i) A proof of (12) can be found, for example, in [2] and in [37]. 
(ii) Using the continuous injection from (G) into the set of continuous, 
bounded functions on G, it is clear that there exists a constant G (depending 

only on G) such that < C||/||hi(g) II/I|J((g) ^ot every / G R^G) and 

all K > 0. Hence, with (5b), it follows that 

lL^(G)n« - C{£{f,f) + ll/llL2(G)n«)ll/llLi(G,i/(dx)) 
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for all / G Dom{£). This is the Nash inequality. It is then possible to apply 
Theorem 3.25 in [8] (see [27] for details) which yields that there exists two 
constants Ki and K2 depending only on n, r and £ such that 

(14) p(t,x,y) < ^ exp{-\ay - ax\ -tK2\a\^), 

where a = (yo ~ XQ)/KtQ for an arbitrary constant i^T > 0, a time to £ (0, 1] 
and fixed points xo,yo ^ [^)^]- For a choice of K large enough in function 
of K2 (hence, K depends only on k and G) and applying (14) with x = xq, 
y = Ho and t = to, we get (13) on the time interval (0,1]. The Chapman- 
Kolmogorov equation can then be used to get (13) on any time interval 
(0,r] for any r>0. □ 

5. Convergence results. Let (a", /o", fe")ngN be a sequence of functions 
satisfying (5a)-(5c). Let (L*^, Dom(L"')) be the differential operator con- 
structed previously, but with (a,p, 6) replaced by (a", /o", 6"). 

Proposition 3. We assume that 

1 L2(G) 1 1 L2(G) 1 , 6" L2(G) b 
^ — , ^ — and ^ — . 



(i) For any a > (and a = if G is bounded), the weak solutions 

of the elliptic PDE [a — L'^)u'^ = f with the Dirichlet (resp. Neumann) b.c. 
converge weakly in H^(G) to the weak solution of (a — L)u = f . The con- 
tinuous version ofu^{x) given by fQga{x,y)f{y)p'^{y)^^dy with g2{x,y) = 
J(^°° e~"^p"'{t,x,y)dt converges uniformly on each compact of G to the con- 
tinuous version of u{x) given by jQga{x,y)f{y)p~^{y)dy, where ga{x,y) = 
J+'^e-^'p{t,x,y)dt. 

(ii) The weak solution of the parabolic PDE J/'^^ = L'^u"'{t,x) with 
the initial condition 99 G L^(G) converges weakly in L^(0, T; Dom(i5)) to 
the weak solution of the parabolic PDE ^"g^''^^ = Lu{t,x) with the initial 
condition (^GL^(G). Moreover, the continuous version ofu'^{t,x) given by 
jQp"'{t,x,y)ip{y)p'^{y)~^ dy converges uniformly on each compact o/M^ x G 
to the continuous version ofu{t,x) given by J(jp{t,x,y)(p{y)p{y)~^ dy. 

Proof, (i) In fact, solving (a — L"-)u"- = / in H^(G, z^(x)) is equivalent 
to solving 



, 2dx\dxJ p"^ J p" 

in H^(G) =H^(G;(ix), with respect to the scalar product of L^(G;dx). Ac- 
cording to Proposition 5 and Theorem 17 in [41], this ensures the conver- 
gence of u" to u in H'^(G). 
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The estimates on p"'{t,x,y) given in Proposition 2 are uniform with re- 
spect to n. This ensures that p'^{t,x,y) converges uniformly on each com- 
pact of ]R!j_ X G^, at least along a subsequence. Indeed, its limit is necessarily 
p{t, X, y). For each compact subset G' of G, it is well known that any weakly 
convergent sequence in H^(G') converges strongly in L'^(G'). Combining all 
these facts allows us to assert that vP converges pointwise to u (see [33] for 
the details). 

(ii) If p" = 1, the weak convergence of in L^(0, T; II^(G)) comes, for 
example, from Theorem 29 in [42], page 101. 

Otherwise, we have first to use the fact that p"(t, x, y) converges pointwise 
to p{t,x,y) in order to deduce that u"-{t,x) = J(jp"'{t,x,y)ip{y)p'^{y)~^ dy 
converges pointwise to u{t,x) = jQp{t,x,y)ip{y)p{y)~^ dy. 

Moreover, it is standard that is uniformly bounded in L^(0, T; II^(G)). 
Thus, converges weakly in L^(0, T; Gj!,(G)) to u. □ 

Let X"' be the process generated by (L", Dom(L")). 

Proposition 4. Set G = R or G = {£,r) and assume that the Neumann 
b.c. are used in the latter case. Let G' = {i' ,r') C G, and r (resp. t^) he the 
first exit time from G' for X (resp. X^). Under the hypotheses of Proposi- 
tion 3, for any starting point x, 

P,o(X",r")-i ^P,o(X,t)"1 

n^oo 

with respect to the topology of C{[0,T];M) x M for all T > 0. 

Proof. The convergence of X^ to X in finite-dimensional distribution 
follows from the convergence of the density transition function in (3), the 
estimates (12) or (13), and the Markov property (see, e.g., the proof of the 
corresponding result in [33]). 

For the tightness of (Pa; o (X„)~^)„gN, according to the Aldous criterion 
[1], it is sufficient to prove that, for any sequence (r"", 5n)n6N such that is 
a .^-'."'-stopping time and (5„ > is deterministic and converges to 0, we have 

p.[|x,v.„-x;^„l>,7]^-^o 

for all ?? > 0. But, with (12) or (13), 

^ymj>v]<[ §exp(-^^)dz, 

JG\iy-v,y+v) € V On J 

where 9 = 1/2 (for G = M) or 6* > 1/2 [for G = {£,r) and Neumann b.c.]. 
Thus, for n large enough so that if/dn > 1, 

sup¥y[\X^n \>v]< Cxb]l'^~^ f exp(-G2z2) dz 
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The tightness of (P^ o (X")^^)„gN follows immediately by application of the 
strong Markov property to Pa:[|X"„^5^ — X"„| > r/]. 

For the convergence of (X",r"), set Q{x) = inf{t > 0\x{t) ^ G'} for any 
continuous function x : ]R_|_ M. Possibly Q{x) = +00 if x{t) stays in G' . It is 

easy to see that © is lower SGiiii-continuous: thcit is, ^ limiiif^ 

if for all T > 0, supig^^r] \x{t) - ^— ^0. 

Let X be a path such that there exists (x"')ngN converging uniformly 
to X on [0,r] for all T > and 0(x) < liminf„_>oo ©(a;")- It is easily seen 
that this means that x remains in the closure of G' between the time Q{x) 
and liminf„^oo 0(2;")- But one knows that this almost surely never hap- 
pens for a trajectory of a one-dimensional regular diffusion process (see [34], 
Lemma V.46.1, page 273, e.g.). Thus, the set of discontinuities of is a set 
of null measure for the distribution P^;. Since r" = Q{X'^) by definition of 
0, it follows that (X"',r'") converges in distribution to (X,t). □ 

6. On diffusions with discontinuous coefficients. In this section we as- 
sume that, if there is a Dirichlet b.c. at I (resp. r), then we extend the 
coefficients over {—oo,£] (resp. [r, +00)) with p = a = l and 6 = 0. This is 
justified by the results of Section 3.1 and Proposition 4. 

If £ > —00 or r < +00 and we are working with Dirichlet b.c, we extend 
the coefficients on M with p = a = l and 6 = 0. 

Let iT" be a (countable) set of points of {i,r), and assume the following 
hypotheses: 

(15a) a,b and p are right-continuous with left-limit, 

(15b) a, b and p belong to C^{[£,r]\J), 

(15c) J is at most countable, 

(15d) there exists e > such that |x — y| > e for any x,y ^ J . 

Let us remark that (15a)-(15d) ensure that the coefficients are of fi- 
nite variation on [l^r\. For a function / satisfying (15a)-(15d), we will 
denote f'{x) the density of the part of its derivative which is absolutely 
continuous with respect to the Lebesgue measure. 

Proposition 5. (i) We assume that G = M. For any given Brownian 
motion B constructed on a probability space (0, (J^t)t>o, J^, P) with {J-t)t>o o-s 
its natural filtration (transformed to satisfy the usual conditions), (L,Dom(L)) 
is the infinitesimal generator of the unique strong solution X to the SDE 



(16) Xt = Xo+ f Jp{Xs)a{Xs)dBs + f v{dx) dLf (X) 
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(17) u{dx) = Y: + (^^^i^M^l + 



with 

'a'{x)p{x) \ dx 

pyxjox -I- y 

and 
(18) 

a[x+) + a[x—) 

(ii) T/ie operator (L,Dom(L)) to^/i t/ie Neumann h.c. at r and i is the 
infinitesimal generator of the unique solution to (16), where v and f3 are 
defined by (17) and (18) and, in addition, I'dr}) = —1 and z^({-^}) = 1. 

Proof. The existence and the uniqueness of a strong solution to (16) 
follow from the results of J.-F. Le Gall in [23] (see also [3]). Thus, it is 
sufficient to prove that the process X generated by (L,Dom(L)) is a weak 
solution to (16). 

(ii) Case of Neumann boundary condition. 

We assume at first that b = 0. 

It follows from the results on Dirichlet forms that a Revuz measure corre- 
sponds to a continuous additive functional under for any xq G [£, r] . For 
x in [i,r], let {Lf {X))t>o be the continuous additive functional associated to 
the Dirac measure 6x at x. We know that, for a Borel measurable bounded 
function g, the process {J^ g{Xs)p{Xs) ds)t>o is the unique continuous addi- 
tive functional corresponding to the Revuz measure g{x) dx. 

Let g be the continuous version of a function g in H^([^,r]). If Id:x i— > x, 
an integration by parts leads to 

rr 

f(Id,5) = ^ a{x)g'{x)dx 

rr 

= \j^ 0''{x)g{x)dx + \Y{a{y+)-a{y-))g{y) 



(19) 



y&J 

+ \a{r)g{r) - \a{l)g{l) 

I 

g{x)nNidx), 



where 



UN = la'{x) dx+'Y ^{a{y+) - a{y-))5y + \a{r)5r - ^a{£)6e. 
y&J 

On the other hand, 

(20) 2£:(Id •/,/)- f (/, Id^) = f a{x)f{x) dx = f f{x)midx). 
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Under F^^ , for any xq £[£,r], the process X can be written Xt = XQ + Alt + -^t > 
where M is a martingale and is a continuous additive functional locally 
of zero quadratic variation. The bracket (M) of M is a continuous addi- 
tive functional characterized by the measure i^m in (20) (see [13], equality 
(3.2.14), page 110), and is 

{M)t= f\{Xs)p{Xs)ds, 
Jo 

so that there exists a Brownian motion B, possibly on an enlarged proba- 
bility space, such that 

ft 



Mt = J^ yJa{Xs)p{Xs)dBs. 



The process N is characterized by the measure pN in (19) (see [13], Theo- 
rem 5.1.3, page 187 and Corollary 5.4.1, page 224), and is then 

Nt = i f\'{Xs)p{Xs) ds + ia(r)Z[(X) - ia(£)L^(X) 



+ Ei(«(y+)-«(y-))^*W- 

Let {Lf{X))t>Q [resp. L^~^{X), L^~(X)] be the symmetric (resp. right, 
left) local time of X at the point x. As both L^[X) and L^[X) are contin- 
uous additive functionals that increase only on {t > 0|Xj = x}, there exists 
a real number ^{x) such that L^{X) =7(x)L^+(X) (see, e.g., [34], Proposi- 
tion 45.10, page 409). 

The occupation time formula for the local time reads 



(21) g{x+)L^+{X)dx = jf g{Xs)d{M)s = a{Xs)p{Xs)g{Xs) ds. 
On the other hand, 

(22) r g{x)Lf{X)dx= t g{Xs)p{Xs)ds. 







As (21) and (22) are true for any measurable and bounded function g, 

1 



a{x-\-) ' 
One knows that 



Lf+(X) - inX) = 2^* dXs. 



Thus, for any x ^ J , 

(23) Lr(X) - Lf-(X) = {a{x+) - a{x-))^{x)L^+ {X) 
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if X G As, by definition, 

Lf(X) = i(Lr(X)+Lr(X)), 

one gets easily tlie value of P{x). 

At point £, Li{X) = -f{£)Ll+{X) and Lf (X) = ^L^+iX). Thus, ^L\X) = 
L^{X). A similar result holds for L[~(X). Hence, X is a weak solution 
to (16). 

If 6/ 0, substituting ae"^ and /oe~*, where ^ is defined in (7), to a and 
p yields immediately the result, 
(i) IfG = R. 

The previous computations can be used with a localization procedure. 
But it is possible to avoid using the theory of Dirichlet forms. With the 
Ito-Tanaka formula (see, e.g., [34], Chapter IV. 45, page 102) and the results 
from [23], X = S~^{Y) is a solution to (16) if Y is the strong solution of 

Yt = S{x) + e~^^''''^^'^^^p{S-^{Ys))/a{S~^{Ys))dBs 
and the proof of Proposition 5 is now complete. □ 

The next results follow from the Ito-Tanaka formula (see [34], Chap- 
ter IV. 45, page 102, e.g.) applied to S{Xt), where S is given by (10), and X 
is the solution to (16). 

Corollary 1. The speed measure m and the scale function S of {X, {¥x)x£g) 
are given in (9) and (10). 

Remark 3. If G = M, this result could also be proved using a smooth 
approximation of the coefficients and the results of [14]. 

7. Approximation by difTusions with piecewise constant coefficients. In 

this section we assume that 6 = 0, and we have seen in Section 2 that it is 
possible to transform a and p in order to remove the drift. Yet, the results 
of this section may easily be extended to the case b^O. 

7.1. The SDEs satisfied by the approximations. To simplify the notation, 
we assume that I > — oo, r = oo and we set, for n G N, J'^ = {xi\£ = Xq < 

< • • • } for some points Xq,Xi, If ^ = — oo and r = +00, we have to 

pick a reference point on M and to use doubly indexed sequences. 

Thus, we set, for / = a and f = p, 

r(x) = 5]i[,.,,.^^)(x)/(xn, 

k>0 

where x^ is a point in [x",x^^;^). 
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Hypothesis 1. For any n G N, the points Xq < < • • • are chosen such 
that J C the minimal distance between two points of ^7" is positive, 
and 



o"||oo + ||p- P"||oo > 0. 

n— >oo 



Since (a,/o) satisfies (15b) and we assume that J C it is clear that 
one may construct, at least on a compact subset of (^,r), such a sequence 
(i7"')neN- For each n E N, the piecewise constant coefficients a" and and 
the set J'^ satisfy (5a)-(5b) and (15a)"(15d), so that the results of the 
previous sections apply. 

Using the occupation time density formula, it follows from Proposition 5 
that the diffusion X" is solution to the SDE 

(24) = + /* ^a^pn{Xn)dB, + J] ~^lLf {X^\ 

fe>0 

where i3 is a Brownian motion, 

~ Q»(x^g+)-a"(xg-) 

a"(x^+)+a"(x^-) ^ ' 

and Lf{X'"') is the symmetric local time of at point x and time t. If the 
Neumann b.c. is used at i, then (3q = 1. If the Dirichlet b.c. is used at i, 
then we consider (24) up to r = inf{t > 0|Xf = £}. 
Let be the piecewise linear function 

^ / _ r) / __71,\ -r) / __rj,\ 



where A;"(x) is such that < x < x^n^^^-jj^i- Then the symmetric Ito- 

Tanaka formula (see [34], Chapter IV. 45, page 102, e.g.; see also [29] for 
a treatment of this case) applied to yields that Y"^ = $"'(X"') is the 
solution to the SDE 

(25) Yr = Yo" + + E f^kLf (Y^, 

k>0 

where = $(x^) and 



./a"(x^+)//9'^(x^+) - Ja''{xl-)/p''{xl-) 

^a^{xl+)/p^{xl+) + ^a^{xl-)/p^{xl-) 

Of course, if the Dirichlet b.c. is used at £, then we consider only Y up 
to r = inf{t > 0|1( = $(^)}. If the Neumann b.c. is used at £, then we set 
(3q =1. The b.c. at r is treated the same way. 
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The infinitesimal generator of is L^" = on (£,+00) \ J'"', whose 
domain Dom(L^") is the closure of the set of continuous, bounded functions 
/ of class on (^,00) \ J'^ such that f{i) = for the Dirichlet b.c. and 
f'{i) = for the Neumann b.c. and 

for each integer k and n. 

Remark 4. li p = 1 and a{x) = a+ on and a(x) = a_ on M_ with 
a+,a_ > 0, then one derives from (25) and (26) that Po[-'^t > 0] = /3 = 
-y/a+/(-y/a+ + -^0^) for any t>0. The geophysical community has already 
noticed that, in a heterogeneous media with a diffusion coefficient taking 
two values a+ and a_ , the parameter P gives the ratio of concentration of a 
fluid in the upper half-space (see [36], e.g.). 

7.2. The skew Brownian motion. Equation (25) shows that the process 
behaves around like a skew Brownian motion of parameter + /2. 
Various points of view and constructions of this process may be found in 
[4, 16, 17, 38]. ... Of course, the skew Brownian motion of parameter 1 (resp. 
0) is the positively (resp. negatively) reflected Brownian motion (this is used 
to deal with the Neumann b.c). 

Let be a skew Brownian motion of parameter 6 G [0, 1] , and set 

(27) r = inf{t>0||Zf| =/)} 

for some p> 0. We will use the following construction of the skew Brownian 
motion, which may be found in [17], Problem 1, page 115: a skew Brownian 
motion can be constructed by flipping the excursions of a reflected Brownian 
motion with a probability 9. 

To be more precise, let i? be a reflected Brownian motion, and {{£"', r'^)}n£N 
be the family of its excursions intervals. These intervals are such that Rgn = 
Rr^ = 0, i?t > on (r,r"), Un6N(^",^") = ^+ and (r,r") n {t,r^) = 
for n^k (see [34], e.g., for the existence of these intervals). Let e"" be the 
excursion attached to the interval n, that is, = R(^t-£")A{r^-i") [note that 
these intervals (^",r"') may not be ordered, which implies that n is just un- 
derstood as a label with a priori no other meaning]. To the excursion n, we 
associate an independent Bernoulli random variable (j„ of parameter 6 with 
value in {—1,1}. The process constructed by 

Zf = cT„er_^„ if tG[r,r"] 

is then the skew Brownian motion of parameter 6. 

The core idea of the algorithm is contained in the following lemma. 
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Lemma 1. The random variables {t,Z^) are independent. Moreover, 
the distribution of r does not depend on (in particular, r is equal in 
distribution to the first exit time of [—p,p\ for a Brownian motion), and 

Fo[z', = p] = e. 

Proof. This is a direct consequence of the previous construction of the 
skew Brownian motion. □ 

As we will see in the sequel, we will need to simulate a realization of Zf 
given {t <t} under Pq for some p> 0, where r is defined by (27). 

Lemma 2. Let R be a reflected Brownian motion. Then, if < i/q < 
yi < P, 

(28) Fo[Z^ £ [yo,yi];t< T]=eFo[Rt g [yo,yi];t< Tp{R)] 
and if -p<yo<yi< 0, 

(29) Fo[Zf €[yo,yi];t<T] = il-e)Fo[Rt€[-yi,-yo];t<Tp{R)i 
where Tp{R) = inf{t > Q\Rt = p). 

In other words, to simulate Z^ , one needs to simulate the position of a 
reflected Brownian motion R at time t given t < Tp{R) [or to reflect the 
position of a Brownian motion at time t given t <t, whose density is given 
by either (39) or (40)], and to use an independent Bernoulli random variable 
for the sign of Zf. 

Proof of Lemma 2. For a trajectory of the skew Brownian motion 
Z^ , set K^", r")}.„gN as the excursions intervals of Z^. For each n S N, set 

= -Z^(^4_£n)/^(j,n_£n) If t € , t"] , which are the excursions of Z^ . We assume 
that <yo < yi < p. Then 

Fo[Zfe[yo,yi];T.pjR)] 



nGP 



\et-in\ G [yo,yi];sgn(e") = 1; 
T_p^p{e"') >t- sup sup |e''(s)|<l 

k s.t. r'=<f" se[0,r'=-£'^] 

where sgn(e") is the sign of the excursion e", and T^p^p{e'^) is the first time 
(possibly infinite) when the excursion e" reaches —p or p. By construction 
(see Section 7.2), the sign sgn(e") of the excursion e" is a Bernoulli random 
variable of parameter 9 which is independent from any other random vari- 
ables involved in this construction. Equality (28) follows easily, and (29) is 
proved in a similar way. □ 
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8. Error estimates. 

8.1. The elliptic case. Set the Dirichlet problem 

(Lu = on[r,£], 
\u{r)=Ur and {£)=U£. 

We study the error of the solution when (a,p,b) is replaced by (a"" , p'^ ,b^) 
in (30) and when 

sup||(a",p",6")(x)-(a,p,&)(x)|| ^0. 



Proposition 6. There exists a constant C depending only on \, A, r 
and (. such that 

sup \u{x) - < C\ur - U(\ sup || (a", p", 6") (x) - {a,p,b){x)\\. 



Proof. We use for L the operator L = 2~^e~*/)^(e*a^). Let a" and p^ 
be some approximations of a = ae* and p = pe~^ , such that a" and con- 
verge uniformly to a and p. 

Let n" be the solution of (30) with L replaced by L", and u the solution 
of (30). 

Then, using f " = — n G Hg(G) as a test function, one gets 

/ a — — dx= I (a"^ -a) — ——dx. 
Jg \ dx J Jg dx ax 

As xy < Xx'^/2 + 2X~^y'^ for all x, y > and all A > 0, it follows classically 
that 



<2A-^||a"-a 



L2(G) 



dx 



L2(G) 



Let if be the linear function if{x) = {u^ — ui){x — t) / {r — tj + U£. Then, u — (p 
belongs to Hq(G), and then 

du\'^ Ur — ue f du . 
I ax = — / a— dx. 



G \dxj 



r — I Jg dx 



Hence, there exists some constant C, depending only on A, r, I, Ur and 
such that 



^Vdx<C 
ax J 



Ur — U£ 
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As n" and u satisfy the same boundary condition, belongs to Hq(G), and 
from the Poincare inequahty, 



\„,n\\2 ^ f-i 
\V IIl2{G) - ^ 



dx 



L2(G) 



Moreover, H^(G) can be continuously injected in the space of continuous 



functions on G. Thus, for some constant C , sup^g^j [^"'(x)^ < C'll-y^Hyi/ 



It follows that 



sup|'u'"(x)-'u(x)| <C7"||a" 
xeG 



Ur — Ui 



r 



for some constant C" depending only on r, i, u^, u^, \ and A. 
To conclude the proof, it remains to see that 

|a"(x) - a(x)| < Ci sup |a"(2;) - a{x)\ 

+ C2 sup I (x) - I + C3 sup 1 6" (x) - I 

xgg xgg 

for some constants Ci, C2 and C3 depending only A, A, r and ^. □ 



8.2. The parabolic case. The parabolic case is harder to deal with, and 
we are not able to give a full treatment of it. 

Let u be the solution of ^ = Lu on M_|_ x G, with the Dirichlet or Neu- 
mann b.c. on M+ X {^,r} and initial condition '^{x). Let also vP" be the 
solution of the similar problem where L is replaced by L". 



Proposition 7. (i) Assume that p = l. Then, for any finite open inter- 
val {i',r') C G, there exists a constant G depending on \, A, \\^\\i,2, £' , r' 
and T such that 

sup |ti"(t,x) <C(||a"-a||oo + ||?)"-?)||oo)- 

(t,x)G[0,T]x(^',r') 

(ii) When p^l, assume that cp belongs to H^(G). Then, for any finite open 
interval {i',r') C G, there exists a constant G depending on X, A, H^jUhi, ^' , 
r' and T such that 

sup |n"(t,x) - u{t, x)\ < C(||/)" - pIIoo + ||a" - a||oo + W - b\\oo). 

{t,x)e[0,T]x(£',r') 



Proof. We set G' 
defined by 



\v\g't - 



',/). For r > 0, we denote by | • \g',t the norm 

1/2 



sup Ib(t,-)||L2(G)' + 

te[o,r] 



\Vv{t, 



Il2(g) 



,dt 
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for any v G C{Q,T;l?{G)') r\l?{Q,T;B^{G)'). The convergence relies on the 
following estimate (which is easily derived from [19] , Chapter II. 3 inequality 
(3.1), page 74 and the Poincare inequality, e.g.): 



for a constant C that depends only on T and G' . 

We now apply (31) to = vP' — u. Indeed, it is easily seen that is the 
weak solution to 



where dtu(t,x) is a priori a distribution in L^{0, T; R~^{G)), where L^{0, T; X) 
is the space of L^-functions with values in a Banach space X. 

(i) In this case, the last term in the right-hand side of (32) vanishes. 
If G is finite, we use u"" as a test function in (32). After integrating with 
respect to t, we see that we are in the same position as in the elliptic case 
and standard computations yield the desired result. 

If G is not finite, we fix £' < r' and we choose a smooth function ^ with 
compact support such that ^(x) = 1 on {£' ,r'). Then, since any function v in 
H^(G) is locally bounded, V(v^) = ^Vv + vVS,. We use v^^ as a test function 
in (32) and the expansion of ^{v^,): this reduces the problem to the case 
where G is finite. Hence, the proposition is proved. 

(ii) As in [19], Theorem III. 6.1, page 178, we show that dtu{t,x) belongs 
indeed to L^(0, T; L^(G)). For that, we assume at first that a, p and h are 
smooth, so that u is also smooth. We transform L into a divergence form 

operator with characteristic (a,p, 0) = (e*a, e~*/9, 0) as in Section 2. We 
assume that G C([(G;M). In this case, using dtu as a test function for the 

PDE dtu = Lu, and integrating by parts against e"^//9 with respect to x and 
with respect to t, one obtains 



(31) 



sup \v{t,x)\ <G\v\g',t 

{t,x)£R+xG' 



(32) 



dv^{t,x) 
di 



LV{t, x) + - V((a"(x) - a{x))Vu{t, x)) 

+ (6"(x) - b{x))Vu{t, x)+(^- -) dtu{t, x) 




(33) 




JG Jg 
As a and p are bounded and positive, (33) yields 



+ / a{x)\ip' {x)\'^ dx - [ a{x)\d^u{T,x)\ 



'dx. 



(34) 
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where C depends only on the constants A and A. 

If the characteristic of L is not smooth and 93 belongs only to H^(G), an 
approximation procedure shows that (34) is still true. 

Hence, with the additional assumption that the initial condition 99 belongs 
to H^(G), it is possible to proceed as in (i). □ 

9. Numerical simulations: the algorithm. In this section we give two 
algorithms to simulate either Xt for a given time or (X^,r) (or Xt given 
by t < t), where r is the first time the process X reaches the points r or £ 
where a Dirichlet b.c. holds. 



9.1. What is computed? . Our algorithm can be used to approximate the 
following quantities: 

• For any initial distribution /i and any t> 0, we may compute the prob- 
ability density function u*{t,y) = jQfi{dx)p{t,x,y) with respect to dy/p{y), 
since u*{t,y) is the solution to the PDE 

^^^-^^ = L*n*(t,y) and u*{t,y) dy-^ fi. 
In this case, we use the approximation, for a given e > 0, 

u*{t,y)^^Card{i = l,...,n\Z'G[y-e,y + e]}, 

where are n independent realizations of Xf for which t> t with the initial 
measure fi. 

• For any fixed x G {£, r) and any fixed t > 0, we may compute the solution 
u{t,x) of the parabolic PDE (11) for any initial condition (p. In this case, 
we use the approximation 

1 

uit,x)^-Y.^iZ'), 

where the 's are defined as previously with the initial distribution 11 = 5x- 
Indeed, it is possible to deal with nonhomogeneous Dirichlet, Neumann 
or Robin b.c. by possibly coupling our algorithm locally with another algo- 
rithm: see Section 9.5. 

• For any fixed x G (£, r), we may compute the solution u{x) to the elliptic 
PDE 

Lu = on (^,r), u{r)=Ur and u{£)=U£ 
for any {uj.,U() G M^. In this case, we use the approximation 
1 " 

k=l 
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where the Y^'s are n independent reaUzations of under P^;. If there is a 
Dirichlet b.c. at one of the endpoints, one may also consider a Neumann or 
Robin condition at the other endpoint, since a probabihstic representation 
of the solution similar to the one given in the parabolic case still holds. 

Besides, if /? = 1 and 6 = 0, then L is self- adjoint with respect to the 
Lebesgue measure, and computing ip{x)p{t, x,y) dx is sufficient to com- 
pute the solution u{t,x) for any x G (£,r), since p{t,x,y) =p(t,y,x) for all 
t > and all x,y £ (r, i) . 

Remark 5. It is possible to get an analytical expression of p{t,x,y) 
when {x,y) belongs to a certain subspace of G (see [15]), but it leads to 
rather complicated expressions involving spectral decompositions. 

If (a,p) are piecewise constant and 6 = 0, then the numerical errors come 
from the following: (1) The approximations of series by keeping a finite 
number of terms (see Section 9.3); (2) The fact that we use a pseudo-random 
generator instead of independent random variables; (3) The Monte Carlo 
error, that is, the replacement of E[y] by N^^ ^k=i where 1" is a random 
variable, and the y*'s are copies of independent random variables with the 
same distribution as Y. 

The errors of type (1) and (2) are very small, and the error of type (3) 
occurs in all Monte Carlo methods and depends on the variance of Y. 

9.2. The algorithm. For a general {a,p,b), we may transform the diffu- 
sion process into a diffusion process with characteristics {a,p,0) and then 
use piecewise constant approximations (a"',p",0) of {a,p,0). An additional 
error comes from these approximations, which is discussed in Section 8. 

Let be the process generated by (a",/3",0). We have seen in Section 7 
that one can find a deterministic bijection such that Yf^ = <I>"(X") is the 
solution to the SDE 



with = $"(x'j!) and /?|? given by (26). 

In this section we give an algorithm to simulate (r A T, ^Tat) for a given 
time r > 0, where r is the first time the process reaches a point at which 
there is a Dirichlet b.c. The computation of (r AT, X^^^") is done by setting 



This algorithm is easily simplified to simulate (Y^,t). 
Around each point y^, y" behaves like a skew Brownian motion of pa- 
rameter /3^. The justification of our algorithm relies on Lemmas 1 and 2. 
For each G <I>"(i7"), one picks two points ?/^'~ and y]^'~^ such that 



Yr = Y,- + Bt + J2Pj:L'[HY-) 




yk-i<yk <yk<yk <yk+i and y^ 



n n n 

yk=yk- yk 



n,— 
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Let us denote by /C" the set /C" = {JkyoiVk'} U Ufc>o{yfe'^}- The sets /C" 
and <I>"(i7") may have commons elements. The points of /C" have been 
introduced in order to use Lemma 1, and then to use random variables 
which are easily simulated. 

The idea is to compute the successive times and positions on /C"" U <!>"( J"") 
of a particle, with a special treatment to deal with the final time. The succes- 
sive positions of the the particle on /C" U $"(^7") give simply a Markov chain 
on a discrete space, but dealing with the time requires a special treatment. 

Let us fix n large enough to get a fine approximation of the diffusion 
process. Since n is fixed once and for all, we omit future reference to this 
integer. 

The following algorithm takes as an input a horizon time T and the start- 
ing point y of the particle. It returns (r AT,Yta,t), where r is the first time 
the particle hits one of the endpoints of [i,r] where the Dirichlet b.c. holds. 

Notation. We denote hy B a Brownian motion and a skew Brow- 
nian motion of parameter f3 £ [—1, 1]. All the random variables simulated in 
this algorithm are assumed to be independent. 

The MAIN LOOP. We now explain how to update the position of the 
particle when it lies at a point y = yk at time t. We start with y = y and 
t = 0, and this loop is executed until the algorithm stops. 

Case y £ ^{J^)- Here, yk corresponds to a point in which the coefficients 
a and p may be discontinuous. 

• If yk is at one of the endpoints of [£, r] where the Dirichlet b.c. holds, then 
return {t,yk) and stop. 

• Let Z^^'^^"^^/^ be a skew Brownian motion of parameter {Pk + l)/2 such 

that = yk. Set r = inf{s > Olzi^'+^^Z^ G {Vk ,yt}} and compute 

7 = Pyj.[T — t <t]. (Note that from Lemma 1, 7 does not depend on Pk 

since yk-Vk =yk - Vk-) 

• Use a Bernoulli random variable of parameter 7 to decide if r > T — f or 
T <T — t (note that r has not yet been simulated) . 

• If T > T — t, then the particle does not exit from [y^, y^] before T: draw 

a realization z of Z^^^^^^"^ given {t > T — t}. Finally, return (T, z). 

• IfT<T-t, draw a realization {t',z) of (r, zl^'+^^Z^) given {r < T - t}. 

Indeed, r and zi^''^^^^'^ are independent, and Zt^'''^^^^'^ is a Bernoulli 
random variable of parameter {/3k + l)/2 with values in {y^,2/^}. Update 
the current position of the particle by setting y ^ z and the current time 
by setting t <^t + t' . Restart at the beginning of the loop. 

Case y ^ ^{iJ)'- Except maybe for the initial point y, this means that 
y £ IC\ ^{J) ■ Indeed, the operations are similar to the previous ones, except 
that one uses a Brownian motion instead of a skew Brownian motion. 
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• Compute -/ = Fy[T-t<T] with r = inf{s > 0\Bs G ^{J)}. 

• Use a Bernoulli random variable of parameter 7 to decide if r < T — t or 

T>T-t. 

• If T > T — t, then draw a realization z of i?T-t given T — t <t under ¥y. 
Then return (T, z) and stop. 

• If r < T — t , then draw a realization z of B-r given t <T — t and afterward a 
realization t' of r given t <T — t and Br = z. Update the current position 
of the particle by setting y ^ z and the current time by setting t ^t + t' . 
Restart at the beginning of the loop. 

9.3. The random variables that should he simulated. Let [B, (Px)xeK) be 
a Brownian motion. Let Ta^b = inf{t > 0|i?r G {a, 6}} for a <b. Using the 
scaling property, o t~1 is equal to IP(x-fe)/(fe-a) ° ((^ ~ ■ Hence, we 

set r = T_i^i and we assume that a = —1, 6=1. 

Mainly, we need to simulate r or some random variables whose distribu- 
tions are related to that of r. Prom a numerical point of view, in order to 
simulate a random variable with distribution function F, we may compute 
F~^{U), where U is a realization of a uniform random variable over [0,1] 
(about the simulation of random variables, see, e.g., [10]). 

We give then explicit expressions of the densities of the random vari- 
ables of interest (see [6] or [28], e.g.). Their distribution functions are easily 
computed, and may then be efficiently inverted using Newton's method. 

If G(t,x) =P^[r <t], then 



(35) 



or 



(36) 



k= 



dG, , 



l + x + Ak _n 4-^4-/1^2/9, l-x + Ak 



(l+x+4fc)V2t _^ ~ "T '^'^ ^-(l-x+4kY/2t 



ttV^, , , /-Tr^(2k + iyt\ f / 1 
= 2 E (- 1) (2^^ + 1) exp ^— '—j cos [xT, [k + - 

Besides, if iJ(t,x) =P^[r <t\Br = l], then 
^^^^ dt^^'^^'-l + x^^ [ V2^t3/2 ^ 

fe = — CXD * 

or 

(38) — (t,x) = V(-l)'=/cexp^ sin — (x + 1) 
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The Bayes formula allows then to compute Px[Bt = 1|t < t] and by symme- 
try, ¥x[Br = -l|r < t]=¥i_x[Br = l|r < t]. 

We also need ¥x[Bt < y\t < r] = F{t,x,y)/¥x[t < r], where F{t,x,y) is 
defined as F(t,x,y) = Fx[Bt < y;t <t]. The density of the Brownian motion 
killed when exiting from [—1,1] is 



or, using a spectral decomposition, 

(40) -Q^it^^^y) = 2 E exp(^ ^ j sm(^ — (x + 1) j sm(^ — (y + 1) j . 

The series giving F, G and H converge very quickly and do not create 
numerical problems. The series (35), (37) and (39) are numerically suitable 
for small time, while (36), (38) and (40) are numerically suitable for large 
time. 

9.4. Efficiency. Most of the computation time is spent in the simulation 
of the exit time from an interval for the Brownian motion. 

If we denote by A^'" = J"^ U /C", where /C" has been introduced at the end 
of Section 9.1, then we need to simulate the exit time of intervals of type 

= K-n-^fc+i] if -^r's are the ordered points of A^". Thus, the "cost" 
of our algorithm may be identified with the number of times one needs to 
simulate the exit time of some I^. 

However, it is difficult to estimate this number N'^ of computations since 
it depends strongly on (a) the distance between z^_-^ and z^_^^ for /c > 1, 
and (b) the number of passages on each interval = for k>l. 

Besides, the number of such intervals depends on the variation of the 
coefficients. Intuitively, the fiatter the coefficients, the more efficient is our 
algorithm. 

In a simple but realistic case, we can give a rough estimate of the cost of 
our algorithm, which is matched by some numerical experiment (see Figure 3 
in Section 10.2). 

We assume that the coefficients a, p and b are of class on M \ {0}. 
Fix n large enough and set A = 1/n. We choose the map <I>" so that the 
set $"(J") = {A;A|fcGZ}. We set (a", p", 6")(x) = (a, p, 6)(fcA+) if xG [A;A, 
{k + 1)A] with A: £ Z. Let u (resp. u^) be the solution of the parabolic PDF 
du/dt = Lu with u{0,x) = f{x) [resp. du'^/dt = L'^vP' with u"'(0,x) = ^{x)], 
where L (resp. L") is a differential operator with characteristics {a,p,b) 




(39) 
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[resp. (o", 6")]. Then, for any {t,x), \u^(t,x) — u(t,x)\ is of order 0(A), 
since the coefficients are Lipschitz continuous on and Ml . Moreover, the 
Monte Carlo error in the evaluation of u"'{t,x) is of order 0{l/y/N), where 
N is the number of random independent particles, which gives a total error 
in the evaluation of u{t,x) of order 0(l/\/iV) + 0{A). 

As the intervals IJ! are of type [{k — 1)A, {k + 1)A], in order to simulate 
Xi, we have to simulate A^"^ independent random variables r^,T^, . . . giving 
the exit time from [—A, A] for the Brownian motion starting from 0, where 

iV^ is such that H h t^^~^ < 1 < H h r^^. Replacing r by its 

average E[r] = A^, one gets that N'^ ~ 1/A^. This means that the cost is of 
order 0(1/A^) for a trajectory. Thus, if one wants \u{t,x) — u'^(t,x)\ to be of 
order 0(A), then one has to choose N = 1/A^. This means that, for a weak 
error of order A, the number of times the random variable r is simulated is 
of order 1/A^ with our algorithm. 

9.5. Coupling with other algorithms. Of course, this algorithm can be 
coupled with other algorithms (such as the Euler scheme, e.g.) if the coeffi- 
cients are smooth enough outside some subset / of {i,r). Let us explain 
our idea: assume that the coefficients are smooth outside some interval 
{i',r'). Then, one can pick some points r" and i" such that r < r" < r' 
and i' < i" < I. One can then use the algorithm which is the most adapted 
to the situation for the simulation X outside {i',r') and this shall be done 
until X hits r' or i' . Inside {i',r') one can use the algorithm proposed here 
for the simulation X until it hits r" or i" . 

One can also use a scheme where the process is killed at the discontinuities 
(there are efficient adaptations of the Euler scheme in this case). Thus, the 
distributions related to the skew Brownian motion may be used to re-inject 
the particle in the media. 

Another way of coupling consists in using locally a scheme to deal with 
non-homogeneous Dirichlet, Neumann or Robin boundary conditions. For 
that, we use the following representations when the lateral function ip : — > 
R is bounded and continuous (to simplify the notation, we assume that 
r = +00 and we set r = mf{t > 0\Xt = i}): 



if au{t, r) + n'(t, r) = ilj{t) and q > 0. 

With a Neumann or Robin b.c. at ^, one has to consider the diffusion pro- 
cess which is reflected at this point. The Lepingle scheme [24, 25] gives an 



u(i, x) 



E,[(^(At);i<r]+E,[V;(r)] 



iiu{t,l) = ^{t) 
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easy way to simulate the couple {\Bt — £\,LI{B)) for any t > 0, where B 
is a Brownian motion. Thus, one may then simulate {Xt+st, Ll_^_g^{X)) by 
this way when Xt = i and 6t is small enough. This gives approximations 
of integrals of type JqiP{s) dLl(X), and allows to compute u{t,x) using the 
previous formulae. 

9.6. Localization. If G = R, it follows from Aronson's estimates (see 
Proposition 2) that 



sup Pa: 



xGlf 



sup \Xs — x\> R 



<Ciexp{-C2R^/t) 



for some constants Ci, C2 depending only on A, A (see, e.g., [37] for a proof). 
Thus, one can assume that the coefficients a, p and b are constant far enough 
from the starting point, or that the process can be killed when it reaches 
the edges of a finite interval rather than dealing with a process that lives on 
the whole space M. 



10. Examples. 

10.1. The doubly skew Brownian motion. We apply our algorithm to 
simulate the density of the solution to the SDE 

Xt = Xo + Bt + Il;'/\x) - lLy\x), 

where B is a Brownian motion. Such a process may be called a doubly skew 
Brownian motion. 

This process is then generated by 

2 dx \ dx J 

with 

(a(x)p(x)) = |(^'^)' if ^^(-1/2, 1/2), 

{a{x),p[^x)) 1(2,1/2), if xG (-1/2,1/2). 

In Figure 1 we represent the density at four different times. As expected, 
the density is more concentrated on the interval (—1/2,1/2), which agrees 
with our intuition, thanks to the choice of the coefficients. 

10.2. Diffusion with a coefficient discontinuous at one point. We con- 
sider now that 6 = 0, p = 1 and that 

, \ J 2 + sin(x) , if X < 0, 

"^^-^ ~ |5 + sin(x + 7r), if x > 0. 

We choose the points of T" such that <I>"(J") = {fcA|A; G Z} with A = 1/n. 
The density p(t,x,y) for t = 1.0 and x = 0.5 of the process is represented 
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in Figure 2 for A = 0.1 and 10.000 particles, while the number of times 
the exit time from an interval for the Brownian motion have been drawn is 
shown in Figure 3. Figure 3 is fully in agreement with the rough estimate of 
Section 9.4. 




Fig. 1. Density p{t,x,y) of the doubly skew Brownian motion with x = 1 and t — 0.5, 
t = 1, t = 1.5 and t = 2 (with 50.000 particles). 




Fig. 2. Densities pit, x,y) of the diffusion of Section 10.2 with t = 1.0 and x = 0.5: The 
full line IS for the histogram obtained with the scheme presented m this article. The dashed 
line IS for the histogram obtained with an Euler scheme. 
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10.2.1. Comparison with the Euler scheme. In Figure 2 we have also 
drawn the density obtained with the Euler scheme presented by one of the 
authors in [27]. Indeed, the process X is solution to the SDE 

JO 2 Jo a(0+) + a(0-) 

Using the function 

/ N 1-2/3 , ^ 1 

with 

/3 = (a(0+)-a(0-))/2a(0+), 

one obtains from Ito-Tanaka that Yf = (j){Xt) is the solution to some SDE 
with coefficients that are discontinuous at 0, but without local time, for 
which an Euler scheme is possible and may be applied. The convergence of 




IJ I I I ! I I I I I I 

Ml 0.02 oxa i.m (j,(k u.«6 u.ot u.ib am u.io o.ii 

Fig. 3. in function of A, where N is the average number of simulations of an 

exit time. 




Fig. 4. Interpolation of the solution u{t,x) of (41) computed at points k/2 with k gZ 
fort = 0.5. 
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Fig. 5. Interpolation of the solution u{t,x) of (41) computed at points k/2 with k £Z 
fort =1.0. 

this scheme with discontinuous coefficients is proved in [39], and [27] provides 
an estimation of the speed of convergence of this scheme in this particular 
case. In Figure 2 we use the time step 6t = 0.01 and we see that the two 
empirical densities agree. 

10.2.2. Comparison with a deterministic scheme. We consider now the 
PDE 

du 

(41) —{t,x)=Lu{t,x) with u{0,x) = (p{x), 

where ^p{x) = cos(x) if |x| < 7r/2 and (p{x) = otherwise. 

With our scheme, we computed u{t, x) at time t = 0.5 and time t = 1.0 and 
at the points k/2 with A; G Z. We also use the one-dimensional solver pdepe 
provided with Matlab to solve (41) (indeed, we use Dirichlet boundary 
conditions at —15 and 15, but this does not affect the computations for 
small times). In Figures 4 and 5, we see that the interpolated curves agree, 
even with a small amount of particles (here, 5.000 were used for each starting 
point). 
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